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Quasiparticles in semiconductors — such as microcavity polaritons — can form condensates in 
which the steady-state density profile is set by the balance of pumping and decay. By taking account 
of the polarization degree of freedom for a polariton condensate, and considering the effects of an 
applied magnetic field, we theoretically discuss the interplay between polarization dynamics, and the 
spatial structure of the pumped decaying condensate. If spatial structure is neglected, this dynamics 
has attractors that are linearly polarized condensates (fixed points), and desynchronized solutions 
(limit cycles), with a range of bistability. Considering spatial fluctuations about the fixed point, 
the collective spin modes can either be diffusive, linearly dispersing, or gapped. Including spatial 
structure, interactions between the spin components can influence the dynamics of vortices; produce 
stable complexes of vortices and rarefaction pulses with both co- and counter-rotating polarizations; 
and increase the range of possible limit cycles for the polarization dynamics, with different attractors 
displaying different spatial structures. 
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I. INTRODUCTION 

The experimental realization of quantum condensation 
of quasiparticles with finite lifetimes has opened the pos- 
sibility to study situations where the steady states of the 
condensate are not just controlled by energetics, but in- 
volve also consideration of pumping and decay, leading 
to steady states with quasiparticle currents. The explo- 
ration of the possible behavior exhibited in these systems 
is driven particularly by recent experimental progress in 
microcavity polaritons 1 -^ 3 -^^ and magnonsi^^ (as 
well as more exotic realizations in helium 11 ). Microcavity 
polaritons are the quasiparticles that result from strong 
coupling between photons confined in planar semiconduc- 
tor microcavities and excitons confined in quantum wells. 
The photons are confined by Bragg reflectors to a two- 
dimensional cavity, and for small in-plane momentum 
have an effectively quadratic dispersion, with a mass of 
the order of 10 -5 - 10 -4 of the free electron mass. From 
strong coupling between these photon states and exciton 
states in quantum wells, new normal modes arise: polari- 
tons. These are hybrid light-matter particles, which have 
both a light effective mass, as well as sufficient interac- 
tions to allow a quasi-equilibrium particle distribution to 
be established. However, the distribution is only quasi- 
equilibrium, since the photon component is not perfectly 
confined, so polaritons have a finite lifetime, and to main- 
tain a steady state population, continuous pumping is re- 
quired. Even when the energy distribution of polaritons 
looks reasonably thermal, this does not mean effects of 
finite lifetime may be neglected, as one may have ther- 
malized energy distributions on top of steady state flows 
driven by the pumping and decay. 

Experimentally, it is observed that above a thresh- 
old pumping strength, an accumulation of low energy 
polaritons is accompanied by: a significant increase in 



temporal coherence, spatial coherence that extends over 
the entire cloud of polaritons^, the appearance of quan- 
tized vortices 3 -, suggestions of changes to the excitation 
spectrum 4 , and a robustness of polariton propagation 
to disorder 5,6 . Of these, the observation of spontaneous 
vortices 3 in polariton condensates provides a strong hint 
that the steady states of the system are affected by the 
flow of particles. For a more general introduction to mi- 
crocavity polaritons, their condensation, and the physical 
systems studied there exist several reviews or books, e.g., 
Refs. HEIII 

The aim of this article is to investigate the rich varieties 
of behavior that come from combining spatial profiles 
driven by particle flow with the dynamics due to the spin 
degree of freedom that microcavity polaritons also pos- 
sess. The effects of this spin degree of freedom have been 
previously considered in equilibrium: due to the weak at- 
tractive interaction between opposite spin polarizations, 
a condensate is expected to be linearly polarized 15,16 , i.e., 
to have an equal density of left- and right-circularly po- 
larized components. Applying a magnetic field can then 
lead to a phase transition, as the density of one circular 
polarization increases, and the other decreases, leading 
through an elliptically polarized phase to a single circu- 
larly polarized condensat o 17 i 18 . In the absence of fur- 
ther symmetry breaking, the linear or elliptically polar- 
ized phase has two gapless modes, corresponding to two 
broken symmetries: the overall polariton phase, and the 
direction of linear polarization (i.e., orientation of ellip- 
tical polarization). The transition to a circular polarized 
state is thus a phase transition, as such a phase has only 
a single broken symmetry, and a single gapless mode. 

Assuming no other symmetry breaking terms, the lin- 
early polarized condensate can have independent vortices 
of the left- and right-circularly polarized components 19 ; 
such polarization-dependent phase winding has recently 
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been observed^. These two vortices are not completely 
independent. Even in the absence of any further symme- 
try breaking, there is a short-range density-density in- 
teraction between vortices of opposite polarizations^*^, 
however such a term does not depend on the circulations 
of the two vortices. If one also takes account of TE-TM 
splitting of the electromagnetic modes 21 , this leads to a 
term that splits linear polarization states according to 
whether the polarization is parallel or perpendicular to 
the polariton momentum; as such, this provides a circu- 
lation dependent interaction between vortices of the two 
polarizations^ 2 -. In addition, real materials are expected 
to possess a small splitting between linear polarization 
states due to asymmetry of the quantum-well interfaces 
in non-centro-symmetric crystals 2 -; such a splitting can 
also be controlled by applying electrostatic fields— or 
stress^ 5 -. Such terms will again induce interactions be- 
tween vortices of the left- and right-circular polarization, 
and if strong enough, will lead to a pinning of the polar- 
ization of a polariton condensate, as has been observed 
in experiment ^ 26 ! 27 ! 28 . 

Considering the splitting of linear polarizations as a 
phase-locking term between the two circular polarization 
components, and combining this with a magnetic field 
that favors one or other circular polarization, one has 
a Josephson problem 29 , where the energy favoring lin- 
ear polarization is a Josephson coupling, and the mag- 
netic field leads to a bias between the two fields. She- 
lykh et al.— have considered the interplay between these 
spinor Josephson oscillations, and Josephson coupling 
between two different spatial modes of a trapped po- 
lariton condensate, showing that complex behavior can 
arise in this four-mode system without pumping and de- 
cay. A similar problem has also been studied in the con- 
text of spinor condensates of cold atoms in double well 
potential o 31 ! 32 . If one does include pumping and decay 
then even the two- mode system (i.e., just the dynam- 
ics of the spin component) can show a rich variety of 
behavio r 33 ! 34 . Due to dissipation, the dynamics settles 
on an attr actor, which may either correspond to a phase- 
locked (i.e., linearly polarized) state, or may correspond 
to a limit cycle, in which the phase difference between the 
two polarization components continually increases, with 
the cyclic nature arising from 2tt periodicity in the phase 
difference. There also exist parameter ranges where there 
is bistability, so that which attractor the system settles on 
depends on the initial conditions. In this article, we will 
both review this two- mode dynamics, and then extend to 
the case of including many spatial modes in addition to 
the spin dynamics. 

To describe multiple spatial modes in a trapped, 
pumped, decaying condensate, a mean-field approach 
to this problem leads to a complex Gross-Pitaevskii 
equatio n 35 ! 36 . It is also possible to include fluctuations 
beyond mean-field theory within such a formalism by 
quantum stochastic approaches 37 . As well as describing 
the spatial structure of nonequilibrium condensates, the 
complex Gross-Pitaevskii equation has also been applied 



in a wide range of areas, see for example Refs. I38ll39l . 
One application which is very closely connected to mi- 
crocavity polaritons is the dynamics of lasers propagat- 
ing through strongly nonlinear materials, where combi- 
nations of cubic and quintic nonlinearities can produce a 
state with a preferred density, sometimes referred to as 
"liquid light" 4Q ! 41 . The results in this article suggest that 
consideration of polarization dynamics in such materials 
may also provide interesting results. 

By considering the spinor complex Gross-Pitaevskii 
equation to describe the spinor condensate in an har- 
monic trap, several new classes of attr actors are seen to 
occur in addition to those present for the two-mode sys- 
tem. There are limit cycles describing small oscillations 
of the phase difference between the two components, and 
limit cycles with 4tt periodicity of the phase difference. In 
addition, the different dynamical attractors correspond 
to different spatial profiles, in which the presence or ab- 
sence of vortex lattices can be influenced by the applied 
magnetic fields. 

As well as describing the steady state profiles, the com- 
plex Gross-Pitaevskii equation can describe the normal 
modes for small fluctuations about such steady states. In 
the absence of a spatial trap, one notable consequence of 
including pumping and decay is that the long- wavelength 
modes of a dissipative condensate are modified, to be- 
come diffusiv e) 42 1 43 ! 44 ; similar results are seen both for in- 
coherent pumped condensates, and for optical parametric 
oscillation, where scattering of pumped polaritons into 
signal and idler states is modified by bosonic enhance- 
ment due to the population of the signa l 45 ! 46 . From this 
standpoint, another purpose of this article is to address 
how the combination of pumping, decay, and symmetry- 
breaking terms modify the spectrum of the spinor con- 
densate. While a diffusive mode for the overall polari- 
ton density and phase always exists, the long- wavelength 
modes for spin modes can show a range of behaviors: dif- 
fusive, linearly dispersing, or gapped, dependent on the 
balance of decay and symmetry-breaking terms. 

The remainder of this article is arranged as follows. In 
Section II we introduce the model of spinor condensates 
as a system of coupled Gross-Pitaevskii equations with 
pumping and decay terms. By neglecting the spatial vari- 
ations of polarized components we discuss the bifurcation 
diagrams for a two-mode system in Section III. The nor- 
mal modes of the homogeneous system are obtained in 
Section IV. In Section V we study the stability of cross- 
polarized vortices. The detailed study of the dynamics 
of the full trapped system is performed in Section VI. In 
Section VII we discuss how various regimes can be de- 
tected in experiments. The conclusions in Section VIII 
summarize our findings. 



II. MODEL FOR SPINOR POLARITONS 

The model we consider in this article consists of the 
complex Gross-Pitaevksii equation (cGPE) as described 
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in Ref. [35|, taking into account the two possible polariza- 
tions of polaritons, written in the basis of left- and right- 
circular polarized states, denoted by i/jl,r- In addition to 
the terms discussed in Ref. 35, three new parameters arise 
from considering the spin degree of freedom in an applied 
magnetic field. Firstly, in addition to an interaction with 
the total polariton density H Uo = (£7 /2)(|V ; l| 2 + \^r\ 2 ) 2 , 
there is an attractive interaction between opposite spin 
species, Hjj 1 = — 2Z7i|^z / | 2 |V ; j r| 2 - Secondly, there is 
a magnetic field term Hq b = (Q B /2)(\^l\ 2 - \^r\ 2 ), 
and finally there is a symmetry breaking term Hj 1 = 
Ji(^l^r + H.c), which may naturally aris o 1 ^ 8 due to 
asymmetry at the quantum- well interfaces 23 , or may be 
induced by electric fields 2 ^ or stress 2 ^. Put together, 
these yield the coupled cGPE: 



ft 2 V 2 
2m 



+ V(r) + U \^l\ 2 + (U - 2/7x)|^| 2 



n L 



+ i (7eff- k-T\i/> l \ 



i^L + Mr, (1) 



with the analogous equation for ijjR following by replac- 
ing ipL <-> i/jr and Qb —> In this article, we 
will either neglect V(r) (in sections IIII| HV]h or con- 



sider a harmonic trap, V(r) 
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case, we will rescale energies and lengths in terms of the 
harmonic-oscillator energy Tiojq, and length I such that 
hujQ = h 2 /2ml 2 . By additionally rescaling the density 
such that U\ip\ 2 — > hcoo\ip\ 2 the cGPE can be written as: 



V 2 + v(r) + |^l| 2 + (1 - u a )\^ R \ 2 + - 



+ i(a-(j\^ L \ 2 ) 



^ L + Ji)R (2) 



with the definitions a — (pf e ^ — n) /Tiuj$,<j = T/Uo, J = 
Ji/huo,A = Q,B/TioJo,u a = 2U\/Uq. If considering a 
harmonic trap, then v(r) = V{lr)/fiujQ = r 2 , otherwise 
we will take v(r) =0. 

The dimensionless parameter u a describes the extent 
of spin anisotropy of the interaction; using the estimate 
in Ref. [47| we take u a = 1.1. From estimates of the polar- 
ization splitting in Refs.[27.28, we take J\ ~ 0.1-0.2meV. 
Using the estimates of trap frequency given in Ref. [35|, 
we will take < a < 10, a = 0.3, and we may write the 
dimensionless splitting J c± 1; in the numerical results, 
we will consider a variety of values of J around this value. 

The spatially extended, harmonically trapped sys- 
tem is unstable unless pumping is restricted to a fi- 
nite spot^. For our studies of this system, we there- 
fore consider a circular pumping spot with a cutoff ra- 
dius ro: a — > aO(ro — r). For constant pumping 
strength, ro determines whether the system forms vor- 
tices spontaneously 3 -^. 



III. REVIEW OF THE TWO-MODE MODEL 

In order to provide a basis from which to understand 
the behavior observed in the trapped system, it is first 
useful to review the results that occur in the "two-mode 
model" , where spatial variation of each component is ne- 
glected, so the coupled complex Gross-Pitaevskii equa- 
tion reduces to two equations for the variables ^l^r- 
This model was studied by Wouters 33 , where the con- 
ditions for the existence of a steady state (i.e., synchro- 
nized) solution were found. 60 As A increased, Wouters 33 
found that a synchronized solution could persist for some 
range of A, with an increasing phase and density differ- 
ence between the two modes, until at some critical A the 
steady state vanished. Near the critical A, the existence 
of bistability was noted. This section will both summa- 
rize these previous results, as well as analyze the behavior 
of the desynchronized solution, showing that the dynam- 
ics is analogous to that of the damped driven pendulum, 
or current-biased Josephson junction 48 . 

To discuss the dynamics, it is convenient to 
reparametrize using: 



^L,R y/PL,Re 



i(<t>±6/2) 



R = 



PL + PR 



PL ~ PR 



and write coupled equations for R, z. (The global phase 
(j) does not affect the dynamics.) In visualizing the dy- 
namics, one may consider a Bloch vector, defined by 
x = ^JpLpR cos(0), y = y/pLpR sin(0), z, in terms of the 
above variables. Since there is pumping and decay, the 
length of the Bloch vector is not conserved, hence there 
is a dynamical equation for R = ^ x 2 + y 2 + z 2 . The 
coupled equations have the form: 



= —A — 2u a z ■ 



2Jz cos(<9) 
VR 2 - z 2 



z = 2(a- 2aR)z - 2J^R 2 - z 2 sin(0) 
R = 2a(^R-R 2 -z 2 Y 



(3) 

(4) 
(5) 



Writing the equations in this form firstly allows one to 
understand the basic effects of pumping and decay on the 
dynamics, and, secondly, will be shown below to reduce 
to the damped driven pendulum in a relevant limit. 

For a steady state, it is clear that the condition R = 
defines a Bloch surface, R(z). One may then study the 
dynamics of points that start on this surface to character- 
ize the nature of the transition as A increases. As shown 
in Fig. [H each starting point z, is attracted either to a 
fixed point, or to a limit cycle. The fixed point is the syn- 
chronized solution; the limit cycle is the desynchronized 
solution, in which the phase difference between the two 
spin components is changing — this behavior is a cycle 
since the dynamics are periodic under 6 — > 6 + 2tt. As A 
increases, there is a lower critical A C) i ower below which 
the all initial conditions on the Bloch surface flow to 
the synchronized solution, and an upper critical A C)Upper 
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FIG. 1: (Color online) Basin of attraction of fixed point and 
limit cycle for full two-mode Josephson problem. Each panel 
is for a given value of A, as indicated, and shows a colormap 
according to the final state found by starting from the given 
initial conditions 6, z. Points that flow to the fixed point are 
colored gray, others are white. The fixed point is marked by a 
red cross, and the limit cycle (when it exists) by the magenta 
line. 



which corresponds to the previously found 33 point be- 
yond which synchronized solutions can no longer exist. 

The behavior seen in Fig.[T]can be understood by notic- 
ing that the choice of parameters used puts one in the 
"Josephson regime" of the Josephson equations^, i.e., 
u a R ^> J, and that in addition the typical dynamics 
obey z <C R. In this case, Eq. ([5]) simply reduces to 
R = a/ a. Then, by eliminating z from Eq. (|3|4p , one 
finds: 

+ 2a0 = -2aA + 4u a J- sin(0). (6) 

a 

This equation describes a damped driven pendulum, or 
alternatively a current-biased Josephson junction. (N.B., 
due to our choice of sign of J, "gravity" for the pen- 
dulum acts to drive 6 ^ tt). The behavior of Eq. (|6]) 
is well known (see, e.g., Ref. |48| ); for large damping 



a there is a simple transition at A c 
between fixed points for A < A 



c, upper 



,er = 2u A J/a 

and limit cy- 
cles above. For smaller damping (explicitly, for for 
a < 0.595 ^/4ii a Ja/<j), there is a range of bist ability, 
where both limit cycles and fixed points may be found^. 
The inset of Fig. [2] shows the bifurcation diagram de- 
termined from Eq. (|3]-[5|), and for comparison, the value 
Ac, upper and the critical damping that result from the 
approximations leading to Eq. (|6|). 

As well as explaining the classes of behavior seen, 
Eq. (|6|) may be used to explain the critical behavior near 



A C) iower- For a sufficiently small that bistability exists, 
then near A C5 i ower , the period of the limit cycle grows 
like T oc | In [A — A Cj i OW er]|- In the context of the cur- 
rent problem, this period relates to the average chem- 
ical potential difference between the two components, 
(0) = (/iL — I^r) — 2tt/T. For large A, one eventu- 
ally reaches \±l — I^r — A. The evolution of (/jLl — I^r) 
is illustrated in Fig. [2j which plots the spectral weight 
\iI>l(w)\ 2 , choosing initial conditions so the limit cycle is 
obtained for all A > A c 
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FIG. 2: (Color online) Fourier transform of ipL(t), showing 
bifurcation of frequencies as A is increased, and illustrating 
critical behavior u ~ 1/ In | A — A c |. Initial conditions for each 
A are chosen so that a limit cycle (i.e., desynchronized solu- 
tion) is found for all A > A C5 i OW er- Inset: Regions of stability 
of the synchronized and desynchronized solutions found from 
Eq. J3H5])- The dotted line marks the approximate A c?up per 
appropriate in the Josephson regime, Eq. (|6]), and the * marks 
the known bifurcation points in that regime. 



IV. NORMAL MODES IN THE EXTENDED 
HOMOGENEOUS SYSTEM 

Before considering the interplay of spin dynamics with 
the spatial profile of a trapped condensate, one may first 
consider a simpler problem involving the interplay of spin 
and spatial dynamics — the finite momentum normal 
modes of the pumped, decaying spinor condensate. The 
normal modes of the spinor condensate without pump- 
ing and decay were discussed i n 17 i 18 . In the equilibrium 
case, for small |A|, there is an elliptical condensate (i.e., 
a finite density of both spin components) and there are 
two gapless linear modes, describing excitations of the 
global phase, and the relative phase of the two compo- 
nents. When |A| becomes large enough to cause a tran- 
sition to a circularly polarized state, only a single gapless 
mode survives, describing phase modes of the condensed 
component. On the other hand, the presence of pump- 
ing and decay is known to replace the linear dispersion 
of the phase modes with a diffusive behavior at small 
momentum 42,43,44,46 . The aim of this section is to see 
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how these two effects are combined. The result is that 
while the global phase mode remains diffusive, the real 
part of the relative phase mode can be either gapped, lin- 
early dispersing, or diffusive according to the value of A 
chosen. In the following, we will first present numerical 
results for the normal-mode frequencies, and then dis- 
cuss how these can be straightforwardly interpreted in 
the same Josephson regime as discussed above. 

For the purpose of numerically calculating the nor- 
mal modes, it is simplest to consider the Bogoliubov 
parametrization of fluctuations at some wavevector k (al- 
lowing for decay) as: 




7i/4 7i/2 371/4 71 7C/4 7i/2 3ti/4 n 




^z,,i?(t,r) = e 
e 



u L,R t 



(7) 



Substituting this into Eq. (J2J), and linearizing in u, v 
yields the secular equation Det[(w — in)l — M] = 0, where 
in the basis \ = (ul, vl, ur, vr) t , the matrix M is 



/ A L 



M : 



-B* L 
A R 
-B*r 



B L 

-A* r 



C L 

-D* T 



-a 



Br Cr Dr 



(8) 



-A* R 



- D* 



Noting that for plane waves, the kinetic energy term 
— V 2 — > /c 2 , the matrix elements are: 



(l-i<7)|<| 2 , B L = {l-ia){^ L f 



C L = J+(l-u a )^° R *, 



D L = (1- UaWl^l, 



and similarly with L <-> R. The normal modes calculated 
this way are shown in Figs. O [H Figure [3] shows the 
modes at k = (lower panels) as well as the value of A, 
and the densities Pl,r (upper panels) corresponding to a 
solution with a given value of 6. [This is to make use of 
the fact that the steady state of Eqs. ([3HS]) can be found 
explicitly for A, Pl,r as a function of (see Ref. 1331).] 

In Fig. [3j the range of 6 is restricted to [0,7r], since 
the range [7r, 2tt] is equivalent to this, after swapping 
L <-» R and A — > —A. Within this range, only values 
> 6 C ~ 7r/2 correspond to stable solutions. As only the 
modes at k = are shown, there is always a zero mode 
corresponding to global phase rotations. The other three 
modes divide into an overdamped mode (largest imagi- 
nary part), and a pair of modes which can be either over- 
damped or underdamped, as seen by the bifurcation of 
the real part at = 1.964. Three values of (correspond- 
ing to three different applied detunings A) are chosen to 
illustrate the overdamped, critically damped, and under- 
damped cases, and the normal modes with finite k are 
shown for these points in Fig. [H In the overdamped case 
(for near 7r/2), both the relative-phase and global-phase 
modes are diffusive at small k. When underdamped, the 
relative-phase mode always has a nonzero real frequency. 
When critically damped, the real part of the relative- 
phase mode has a linear dispersion. We will next discuss 



FIG. 3: (Color online) Steady states, and damping of uniform 
fluctuations. All panels are plotted as a function of the phase 
difference, 0, between the two circular polarization compo- 
nents. Top left: densities of components in steady state so- 
lution. Top right: Detuning corresponding to given solution. 
Bottom: Real (left) and imaginary (right) parts of frequency 
for small fluctuations about the steady state. The gray shaded 
region is unstable to small fluctuations. The crosses mark the 
conditions used for the finite k spectra in Fig. [3] Plotted for 
a 6.5, a 0.3, J = 1.0, u a 1.1. 




FIG. 4: (Color online) Spectra of normal modes in a uni- 
form system; showing real (left) and imaginary (right) parts 
of each mode separately. The three rows correspond to the 
three points marked by crosses in Fig. [3] In all cases the to- 
tal density modes are diffusive, while the spin mode changes 
from diffusive (top row, nearest 7r/2), to linearly dispersive 
(middle row) to underdamped spin oscillations (bottom row, 
nearest tt). 



how this behavior can be understood in the Josephson 
regime. 

To understand the modes qualitatively, it is clearest to 
work in the basis of i?, z, 9, (j). Extending Eqs. ([3HS]) to 
include spatial gradients via a Madelung transformation, 
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one finds the equations: 

R=-V (2RV(j) + zVO) + 2aR - 2a(R 2 + z 2 ) (9) 
z = - V (2zV0 + i?V0) + 2(a - 2aR)z 

-2Jy / R 2 -z 2 sin (9 (10) 

• = X^+X* _ (2 _ Ua)R _ JR cQs{6) (n) 



V# 2 - ^ 2 

2Jz 



VR 2 - z 2 



cos(<9) - A, (12) 



where we have introduced the shorthand Xl r = 
(V 2 v /CT/^p^ - (W L ,^) 2 , with L>jR = ± 0/2. 
When considering fluctuations about the uniform state, 
the gradient terms simplify, as expressions such as (V#) 2 
are necessarily second order in the fluctuations, and 
so may be neglected. We will make the same ap- 
proximations as led to Eq. (|6|) (i.e., assume J <C 
u a R,z <C R), and use the steady-state solution in this 
limit (R = a/a, J sin(0) = — az) to simplify the equa- 
tions for fluctuations. If one writes the fluctuations as 
X = (SR, 5z, (50, 56) T then since these are real variables, 
the normal modes have the time dependence = 
Xexp(— icot — hit), giving a secular equation: 

Det [{iuj + k)1 + M + /c 2 Mi + C>(/c 4 )] = 0, (13) 

where M n indicates an expansion in powers of momen- 
tum, to understand the small k dispersion. The first two 
terms are: 



/ 



Mo 



Mi 



-2a 
-2az 
-2 + u a 







-1/2R z/2R 2 
z/R 2 -1/R 




(14) 



(15) 



The nature of the k = normal modes is immedi- 
ately clear: <p has no restoring force, and so has a zero 
frequency oscillation; R has a damping rate 2a, and so 
describes a decaying mode; z and are mixed, and have 
modes with frequencies given by: 

{iuo + K)(iv + K- 2a) - 4u a JRcos(6) = 0. (16) 

Noting that cos(#) is negative, and that the prefactor is 
the same as that of sin(0) in Eq. (|6]), one may write: 



4u a JRcos(0) 



Q 2 



(17) 



where fl p (0) is a "plasma frequency" describing the 
restoring force as a function of angle. The normal modes 



are uj — in 



-ia ± — a 2 . The transition between 



under- and overdamping occurs because for ~ tt the 
restoring force is large so Vt p > a, while as — > 6 C c± n/2, 



the restoring force vanishes, and so there is an interme- 
diate value of (and hence A) at which Q p = a, de- 
scribing the critical damping. Note that whatever the 
choice of parameters, an over damped regime will always 
exist since Q p must always tend to zero at C . However 
if a is sufficiently large the underdamped regime may 
vanish. The approximate expressions for the eigenvalues 
qj — ik = 0, — ia ± ^ft p (0) 2 — a 2 , —2ia match the general 
form seen in the bottom panel of Fig. 03 it is apparent 
that the full problem has some additional variation of 
damping rate with 0, not captured by the approxima- 
tions used in the above expressions. 

One may now explain the linear dispersion at this criti- 
cal damping as arising from degenerate perturbation the- 
ory, allowing a k 2 perturbation to give a oc k splitting. 
Setting Q p = a, and writing uj — ik = —ia + 77, then 
expanding Eq. ([T3]) to leading order in 77, k gives: 



= r] 2 



2u a R 



2u a R 



(18) 



hence describing modes uj — in = —ia + c e ^k. Note that 
although these modes have a linear dispersion of the real 
part as a function of /c, they have a lifetime that re- 
mains finite as k — > 0. As such, this may provide an 
example where linear dispersion of a given mode in a 
condensed system need not imply superfluidity of the 
associated density. Further work is required to deter- 
mine the current-current response function in this sys- 
tem, which will determine whether there is a difference 
of transverse and longitudinal response for spin currents, 
however superfluidity of the spin current is not expected 
here, since the spin orientation is locked by the Joseph- 
son term J. The linear dispersion arises only when the 
damping matches frequency of this phase-locking term, 
producing critical damping. 



V. STABILITY OF CROSS-POLARIZED 
VORTICES WITH PUMPING AND DECAY 

When an harmonic trap is introduced, the existence of 
steady state currents can destabilize the Thomas-Fermi 
like density profile, and lead to spontaneous rotating vor- 
tex lattices^. Before addressing how this instability in- 
teracts with the polarization degree of freedom, we first 
discuss the simpler question of the dynamics and stabil- 
ity of individual vortices of opposite polarization. One 
should note first that a single vortex in a pumped decay- 
ing condensate already has a more complicated structure 
than the vortex without pumping and decay — the vortex 
becomes a spiral vortex, with both radial and azimuthal 
currents; see appendix [K\ for further discussion. 

Considering a spinor condensate system with vortices 
in both polarization components, there is an interplay 
between the weak attractive coupling that invites the 
vortices to have their cores aligned irrespective of their 
circulation, the Josephson coupling that discourages the 
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alignment of vortices of opposite circulation and the de- 
tuning that stabilizes such alignment. To illustrate these 
possibilities we consider several examples. The notation 
(±n, ±m) refers to condensates with vortices of topolog- 
ical charge ±n (±m) in left (right) condensate. From 
Eq. ([2]), with v(r) = r 2 , the following set of stability 
scenarios are found: 

J = 0: All (±1,0) and (±1, ±1) vortex complexes are dy- 
namically stable. 

J ^ 0, A = 0: Solutions (±1, +1) are stable, (±1, 0) and 
(+1,-1) are unstable. Depending on the strength 
of J, vortices may start precessing around the cen- 
ter of the trap, move beyond the boundary forming 
either a pair of rarefaction waves with opposite ve- 
locities or a complex (±1, +1), or disappear at the 
condensate's edge and bringing in two aligned pairs 
of vortices (±1, +1) and (—1, —1). Fig. [5] illustrates 
these possibilities for a = 4.4, ro = 3, a = 0.3 and 
J = 0.5, 1, 1.5, 2. 

J / 0,A ^ 0: For a given J, any sufficiently large A al- 
lows the vortex complexes (±1, —1) and (±1,0) to 
stabilize. Fig. [6] shows such stabilized (±1,-1) 
complex for J = 1, A = 8. 

The stationary state shown in Fig. [5] for J = 1 is 
the vorticity-free, but not radially symmetric, solution 
of Eq. ([2]). These are analogous to the rarefaction soli- 
tary waves of the nonlinear Schrodinger equation discov- 
ered by Jones and Roberts^. In trapped 2D condensates 
these waves were also found 51 : they form as two vortices 
of opposite circulation disappear at the boundary of the 
condensate. In the conservative GPE these waves propa- 
gate with velocities exceeding the velocity of any vortex 
pair. In spinor damped/driven condensates two such so- 
lutions induce flow in opposite directions forming a sta- 
tionary complex. The plots of the real and imaginary 
parts of and i/jr shown in the left panel of Fig.[7J and 
the density plots of Fig. [5] can be compared to the bottom 
panel of Fig. 9 and the left panel of Fig. 10 of Ref. [5l|. It 
also follows from simulations that ipL(x,y) = iPr(x, —y), 
so the found stationary state satisfies a one component 
Ginzburg-Landau equation 

idtij) = [-V 2 ±r 2 ±|Vf^ 

±i(aB(r - r) - a\^\ 2 )]^ + Ji/>(x, -y). (19) 

This also suggests a way to generate solitary waves in 
spinor condensates. The dark soliton (obtained by phase 
imprinting, for instance) will undergo a transverse snake 
instability and form two stationary pairs of vortices of 
opposite circulation, whereas starting with a (±1,-1) 
complex one will obtain stationary rarefaction pulses. 

As A increases, the rarefaction waves in two compo- 
nents lose their antisymmetry and the density minima 
move away from the center as the right panel of Fig. [71 
illustrates. For intermediate A, complexes combining 
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FIG. 5: (Color online) The outcome of instability of the vor- 
tex state (±1,-1) for a = 4.4, a = 0.3, r = 3, A = 0, 
and J — 0.5 (top row), J — 1 (second row), J — 1.5 (left 
bottom) and J — 2 (right bottom). The initial state is 
ip± = \/®([tTF - r 2 )(x ± \y)/yjr 2 + 1/(/xtf£) where ip R = 
ip-,ipL = ^+5 Mtf = 3a /2a and S is the vortex core param- 
eter; see Appendix [A] Depicted are the density plots and 
streamlines of i^l (left panels in first and second row, bottom 
row) and ipR (right panels in first and second row) . Luminos- 
ity is proportional to the magnitude of the density for density 
plots and to the velocity for streamlines. The size of the 
pumping is shown as a circle of radius ro = 3. For J = 0.5 
the final state consists of two precessing vortices of opposite 
circulation. For J — 1 the system evolves into a pair of two 
gray solitons. For J — 1.5 the vortex of the negative circula- 
tion in ipR leaves and a vortex of positive circulation enters 
the cloud. For J = 2 both vortices leave the condensate and 
two pairs of vortices of opposite circulation enter. For J = 1.5 
and J = 2 ipL = ipR, so only ipL is shown. 



a rarefaction wave in one component with a precessing 
vortex in the other component emerge; for even larger 
A two vortices of opposite circulation move around in 
the condensate. Finally, for even larger A, their cores 
overlap and vortices move to the center of the conden- 
sate. These possibilities are illustrated on Figs. [6] and [51 
The cross-polarized vortices deform the condensate to a 
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FIG. 6: (Color online) Stable state of cross-polarized vortices 
(+1, -1) for a 4.4, a 0.3, r = 3, A = 8, and J = 1. De- 
picted are the density plots and streamlines of ipL (left) and 
ipR (right). Luminosity is proportional to the magnitude of 
the density for density plots, and velocity is shown by stream- 
lines. The size of the pumping is shown as a circle of radius 
ro — 3. 
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FIG. 7: (Color online) Real (thick lines) and imaginary (thin 
lines) parts of ipL (solid lines) and ipR (dashed lines) of the 
rarefaction wave complex as a stationary solution of Eq. J2} 
with v(r) r 2 for a 4.4, a 0.3, J = 1 and A = (left) 
and A = 0.4 (right) plotted along the axis of density sym- 
metry. The functions are antisymmetric with respect to the 
origin ^l(0, —y) — iPr(0, y) for A = 0, but not for A / 0. 



slightly oblate form as seen in the density plots of Fig. [6j 
The streamlines and the density plots indicate that the 
vortices coexist with rarefaction waves that rotate in the 
counterclockwise direction. 

The vortex trajectories with pumping, decay and a spin 
degree of freedom are nontrivial. For a one-component 
conservative GPE, a single vortex moves perpendicularly 
to the background density gradient due to the Magnus 
forc o 52 i 53 (the speed of the motion is however nonuni- 
versal, and depends on the global condensate shaped). 
For a two-component conservative GPE with vortices 
in both components, there is an additional advection 
of each vortex by the flow pattern of the other compo- 
nent. Including also pumping and decay, the trajectories 
of the vortices are yet more complicated, as illustrated 
in Fig. [9] for a (+1,-1) complex with J = 1,A = 4. 
Both vortices move along trajectories closely resembling 
epitrochoids, and the distance between the vortices varies 
quasi-periodically with time. Similarly complicated cy- 
cloid trajectories of vortices are known for two-layer fluids 
with one vortex in each layer — such behavior has been 
seen for example in models of tropical vortices 54 . 
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FIG. 8: (Color online) Density plots \iPl\ 2 (top row) and 
I^kI 2 (bottom row) for the outcome of instability of cross- 
polarized vortices obtained by numerical integration of Eq. [2] 
with v(r) — r 2 , a = 4.4, a = 0.3, J = 1 and various A. The 
initial state is the same as in Fig. [5] Luminosity is propor- 
tional to density. The size of the pumping is shown as a circle 
of radius ro = 3. For A = 0.4 two rarefaction pulses form 
a stationary complex. For A = 2 a vortex of negative circu- 
lation in the R component is coupled to a rarefaction pulse 
in the L component. For A = 4 the two vortices of opposite 
circulation do not align their cores. For A = 2 and A = 4 
the complexes precess around the center with the individual 
vortices moving along an epitrochoid; see also Fig. [9] 
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FIG. 9: (Color online) Trajectories of vortices in the com- 
plex (+1,-1) obtained by numerical integration of Eq.[2]with 
v(r) = r 2 , a = 4.4, cr - 0.3, J - 1 and A = 4, and with 
pumping in a radius ro = 3. These parameters are the same 
as for the density profiles in the right hand column of Fig. [8] 
The trajectory of the vortex of positive (negative) circulation 
in the left (right) component is shown as thick black (thin 
blue) line. The start (end) points of the time interval plotted 
are shown as red (green) filled circles. Note that the vortex 
trajectories are not closed. 



VI. THE TRAPPED SYSTEM 

We study the time-dependent problem defined by 
Eq. ([2]) with a harmonic trapping potential v(r) — r 2 
by numerical integration, using a fourth-order, finite- 
difference approximation in space and a fourth-order 
Runge-Kutta method in time. Pumping is restricted 
to a circular spot of radius ro centered at the bottom 
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of the trap, as described at the end of section [III In 
all of the following, we take the experimentally realis- 
tic parameters^ 5 -* 5 - 5 - a = 4.4 and a = 0.3. With this 
choice of parameters, the critical radius can be estimated 
from the Thomas-Fermi approximation for the problem 
without polarization 3 - 5 - as ro c — 4.7. This section will 
present and discuss the transition between synchronized 
and desynchronized states, and the interplay with the 
instability to vortex lattice formation. Section lYl Al first 
addresses the value of A for which a transition occurs, as 
one changes the size of the pumping spot, and coupling 
J. Section IVlBI will then discuss in greater detail the na- 
ture of the new attractors near the critical A that occur 
when one includes spatial degrees of freedom. 



A. Critical A, bistability and internal Josephson 
effect 

Let us start with the simplest trapped problem, taking 
J = 1.0 and ro = 3.0. We study the behavior of the sys- 
tem as the detuning A is increased in the following way: 
starting at A = from a Gaussian initial state, we let 
the system evolve until a steady solution is reached. We 
then increase A in steps, taking for each new A the final 
state at the previous A as the new initial state. For each 
value of A, the chemical potentials of the two polarization 
components are calculated (see appendix ICl for details). 
The result is shown in the left panel of Fig. [TO] The com- 
ponents stay synchronized, sharing a common chemical 
potential, up to A ^ 7.0. As the detuning is increased 
further, the components desynchronize and AC Joseph- 
son oscillations occur between them (Fig. ITT]) . Note that 
as the components desynchronize, the chemical poten- 
tials show an oscillatory time dependence. Fig. ITOl shows 
a time average in this case. 
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FIG. 10: (Color online) Chemical potentials of the two polar- 
ization components plotted against stepwise increased detun- 
ing A. For solutions with time-dependent chemical potentials, 
a time average is shown. 

In terms of a suitably defined, spatially averaged phase 
difference between the components (see appendix |C]) , 
the desynchronization transition corresponds to a transi- 
tion from a fixed point to a limit cycle in the (#, 0) plane, 
in direct analogy to the two-mode problem (see Fig. [12] 
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FIG. 11: (Color online) Josephson oscillations in the desyn- 
chronized solution at A = 12.0, ro = 3.0, J — 1.0 in the limit 
of large time. Left: n = f \ijj\ 2 d 2 r as a function of time. 
Right: The chemical potentials of the two polarization com- 
ponents. (Time given in units of 2/ujq. Blue and red color 
denote L and R components, respectively.) 
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FIG. 12: (Color online) Phase portraits for the trapped prob- 
lem with J — 1.0 and ro = 3.0. Left: synchronized solution 
for small detuning. Right: desynchronized solution in the 
limit of large detuning (corrected for numerical precession). 

and compare with Fig. [1]). 

The right panel of Fig. [10] shows the chemical poten- 
tials of the two polarization components as the detuning 
is increased stepwise in a system pumped in a spot with 
radius ro = 4.0. The progression from the synchronized 
solution at small detunings to the desynchronized solu- 
tion in the limit of large detunings is now more com- 
plicated than for a small pumping spot, ro = 4.0 is 
too small a pumping radius for vortices to form spon- 
taneously in the corresponding one-component system 35 . 
Accordingly, the spinor system with zero detuning de- 
velops circularly symmetric densities, identical in both 
components (but with different phases for the two com- 
ponents). For small enough detunings, the system stays 
synchronized, the densities adjusting to accommodate 
the common, constant chemical potential (left-most pan- 
els of Fig. [13]). 

As the detuning increases beyond A ~ 3.0 the rota- 
tional instability of the high-angular-momentum modes 3 - 5 , 
reappears, and a four- vortex lattice forms in both compo- 
nents. At this point, the smaller component is sufficiently 
depressed by the detuning that the effective critical ra- 
dius is smaller than the pumping radius. The suppressed 
component becomes rotationally unstable and drags the 
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FIG. 14: (Color online) Density profiles along the x-axis 
for representative examples of synchronized and desynchro- 
nized solutions. Left: synchronized solution without vortices 
(ro = 3.0, A = 4.0). Middle: synchronized solution with 
vortex lattice (ro = 6.0, A = 4.0). Right: desynchronized 
solution without vortices (ro = 3.0, A = 12.0). Dashed and 
solid lines indicate left and right-polarized components, re- 
spectively. J = 1.0 for all three examples. 



FIG. 13: (Color online) Progression from fully synchronized 
solution via formation of vortex lattice to fully desynchronized 
solution at ro = 4.0 and J — 1.0. Top: phase portraits for the 
phase difference (a numerical precession has been removed 
from the A = 20.0 plot). Middle and bottom: total densities 
of the L and R components, respectively. 



larger component along. This transition to a vortex- 
lattice solution causes the (common) chemical potential 
to drop, as shown in Fig. [TOl 

While the two components remain on average phase 
locked, the drop in chemical potential at the formation of 
the vortex lattice is accompanied by a time-dependence of 
the chemical potentials, which oscillate together around 
a common mean. Due to the amplitude of the oscillations 
being slightly larger in the larger component, the fixed 
point in the (0, 6) plane turns into a small limit cycle, 
but with only small variations of unlike the 2tt periodic 
cycles in the desynchronized phase. (Middle panels of 

Fig.ua) 

Finally, as A becomes large, the components desyn- 
chronize completely. This situation is shown in the right- 
hand panels of Fig. [131 The difference in space- averaged 
phase now traces out a limit cycle winding through the 
full 2tt. 

The overall progression from synchronized solutions at 
small A to desynchronized solutions with Josephson os- 
cillations holds very generally for different values of the 
Josephson coupling J and also for both small (no vor- 
tices) and large (vortex lattice) pumping-spot sizes. The 
mechanism in both cases is the same: the synchronized 
solution is upheld by adjusting the densities and a steady 
interconversion current forms. After desynchronization, 
the densities revert to profiles as for the single compo- 
nent condensate, with the average density set by the bal- 
ance of pumping and decay, but with an additional time- 
dependent interconversion current. Fig. [14l shows density 
profiles of typical synchronized solutions with and with- 
out vortices, and also a snapshot of a desynchronized 
solution. 

We also performed calculations where we reset the ini- 



tial conditions to a Gaussian density profile, with equal 
phase of the two polarization components for each value 
of A. Based on the results of the two- mode model shown 
in Fig. [TJ such initial conditions would be expected to 
find the limit cycle whenever it exists, whereas the step- 
wise increase of A discussed above is intended to follow 
the fixed point. In these calculations with resetting of 
initial conditions, desynchronized solutions generally de- 
velop at smaller A, as one should expect if there is a re- 
gion of coexistence of synchronized and desynchronized 
solutions. Fig. [15] shows the highest A yielding synchro- 
nized solutions for different ro and different J, both for 
stepwise increased A (top left) and calculation directly 
from Gaussian initial conditions (top right). Comparing 
the two plots gives an estimate of the region of coexis- 
tence. An example (J = 1.0, ro = 3.0) showing the re- 
gion of bistability is given in the bottom left panel. Note 
that the difference between stepwise increased A and re- 
set to Gaussian initial condition for each A vanishes for 
large detunings. This need not be the case for larger ro 
or larger J, where more than one desynchronized, meta- 
stable state (distinguished by numbers of vortices) may 
be possible at high A. 

Our calculations suggest that the spatially extended 
system allows for a rich variation of behaviors. Notably, 
interconversion due to the Josephson coupling between 
the components need not be uniform in space. The 
interconversion rate is easily calculated from the time- 
dependent local phase difference 0(r) = ^r(t) — 0z,(r) 
as 



dp L (r) 



dt 



dp R (r) 



dt 



= 2Jy/p L (r)p R (r) sin(0(r)), 

(20) 

familiar from the Josephson effect as found in any 
textbook^^. Note that Eq. ([4]) describes the same 
physics for the two- mode model. Fig. [TBI shows a map of 
the interconversion rate dtpi\j at two close points in time 
in the large-time limit of the solution directly from Gaus- 
sian initial conditions with J = 1.0, ro = 6.0, A = 16.0. 
This solution is desynchronized and exhibits Josephson 
oscillations. However, at any given time there is inter- 
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FIG. 15: (Color online) Top: Highest A yielding synchronized 
solutions plotted as a function of ro for J — 0.5, 1.0 and 2.0. 
Left: A increased stepwise. Right: resetting initial conditions 
for each A. Bottom left: bifurcation of chemical potentials 
(plotted as in Fig. I10|) for the two approaches, indicating the 
region of bistability for the example J — 1.0, ro = 3.0. Bot- 
tom right: sketch of regions of stability of synchronized and 
desynchronized solutions for ro = 3.0. Compare with inset of 

Fig. m 
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FIG. 16: (Color online) Spatially nonuniform interconversion 
rate dtpi\j at two different, but close, times. J = 1.0, ro = 
6.0, A = 16.0. 



conversion in both directions at different points in space. 



Calculations directly from Gaussian initial conditions 
at large detuning with a large pumping spot (ro = 6.0) 
indicate that solutions are possible in which the polar- 
ization components develop counter-rotating vortex lat- 
tices. This leads to rapid density modulations, particu- 
larly around the edge of the cloud, and a corresponding 
pattern of interconversion in opposite directions. In this 
case, the Josephson oscillations are suppressed and the 
total (scaled) mass n = J \?p\ 2 d 2 r exhibits rapid, small- 
amplitude oscillations. These effects are shown in Fig.fTTl 
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FIG. 17: (Color online) Top from left to right: scaled mass 
n — f \ip\ 2 d 2 r, chemical potentials and interconversion rate 
dtpi\j for a desynchronized solution with counter-rotating 
components. Bottom: densities of the L and R polariza- 
tion components, respectively. (Note that the gray scale has 
been adjusted for each component separately to accentuate 
the density modulations.) J = 1.0, ro = 6.0, A = 32.0. 



B. Phase portraits of more complicated attractors 



The spatial degree of freedom allows a number of be- 
haviors of the (0, 9) phase portrait that are not possible 
in the two-mode problem, which only allows fixed points 
(synchronized solutions) and limit cycles with winding 
number 1 (desynchronized solutions). In the spatially ex- 
tended system, these behaviors are exemplified in Fig. [12] 
for a system without vortices. Both these two classes of 
attr actor can also be seen when a vortex lattice exists. 
In addition, the spatial degree of freedom gives rise to 
several new classes. 

We have already noted in Fig. [13] an example of a syn- 
chronized limit cycle. This can be distinguished from 
the desynchronized limit cycle by the winding number, 
J 6dt/2it, over one period. The synchronized limit cy- 
cle has winding number 0, and the desynchronized cy- 
cles have winding number 1. We also find an example 
of limit cycle with winding number 2, shown in Fig. [18] 
This solution also exhibits another behavior not possible 
in the two-mode model: a retrograde loop. This solu- 
tion appears in calculation directly from a Gaussian ini- 
tial condition at A = 6.4, which is barely above A c for 
J = 1.0, ro = 3.0. For larger A the loop quickly becomes 
a cusp and then disappears. 

The bottom panels of Fig. [18] show two other behaviors 
found when stepwise increasing the detuning in a system 
with a strong Josephson coupling (J = 2.0). Both these 
solutions are basically synchronized (the time average of 
the chemical potential is the same for both components). 
However, the large detuning causes the chemical poten- 
tials to differ at most instances in time, resulting in be- 
haviors similar to the limit cycles with winding number 
0, but which appear to have a chaotic attr actor and/or 
quasi-periodic behavior. 
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FIG. 18: (Color online) Top: Limit cycle with winding num- 
ber 2 and retrograde loop. Left panel shows a blow-up of 
the loop. A numerical precession has been removed from 
these plots. (J = 1.0, r = 3.0, A = 6.4.) Bottom left: 
Possibly quasi-periodic behavior similar to a limit cycle with 
winding number ( J = 2.0, ro = 3.0, A = 14.0 from step- 
wise increase). Bottom right: Possible chaotic attractor 
(J = 2.0, ro = 3.0, A = 10.0 from stepwise increase). 



VII. EXPERIMENTAL SIGNATURES 

To directly observe rotating vortex lattices in exper- 
iments would require time-resolved measurements on 
timescales of the order of the trap frequency, which would 
be challenging with current experimental configurations. 
It is however possible to see signatures of a vortex lat- 
tice in the momentum- and energy-resolved photolumi- 
nescence spectrum, which can be directly measured in 
the far field. The spectral weight is given by the modu- 
lus squared of the Fourier transform of the wavefunction: 



I(tu,k) 



d re 



-ikr 



dte- iut tp(v,t) 



(21) 



As an illustration of this, Fig. [19] shows the spectral 
weight as a function of (uj,k x ,k y = 0). The vortex lat- 
tice as well as the desynchronization transition can be 
seen in the dispersion spectrum. Fig. [19] shows the spec- 
tra of the two polarization components for two different 
solutions. The top panels show a synchronized solution 
that exhibits a vortex lattice without a central vortex. 
Each ring of vortices in the lattices shows up as a side 
band in the spectrum. If there were no vortices, only 
the bottom band would be present. As a contrast, the 
bottom panels show a desynchronized solution that has a 
vortex lattice with a central vortex. The presence of the 
central vortex means that the spectral weight vanishes 



at k x = k y = 0. Note also that the desynchronization 
causes the spectra of the two components to be shifted 
relative to each other, whereas in the synchronized case, 
the two spectra are identical. 
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FIG. 19: Spectral weight at k y = for two vortex- lattice so- 
lutions. Left and right panels show left and right polarized 
components, respectively. Top: synchronized solution with- 
out a central vortex. (J = 0.5, A = 3.5, ro = 6.0.) Bottom: 
desynchronized solution with a central vortex. ( J = 0.5, A = 
16.0, r =6.0.) 



VIII. CONCLUSIONS 

The interplay between the effects of pumping and de- 
cay in setting the density profile, and the dynamics of 
the polarization degree of freedom lead to a rich vari- 
ety of possible nonequilibrium steady states, as well as 
dynamical attractors for the polarized polariton conden- 
sate in a harmonic trap. Even neglecting spatial currents, 
the two-mode model shows nontrivial behavior as a func- 
tion of applied magnetic field: with strong damping there 
is a simple transition between synchronized and desyn- 
chronized states, whereas for weak damping, a region of 
bistability exists, in which the state established depends 
on the initial conditions. 

Allowing for spatial fluctuations in a homogeneous sys- 
tem, the plane- wave modes on top of the synchronized so- 
lution show a second class of transition: if the phase lock- 
ing term for the spin degree of freedom is large enough, 
then at weak magnetic fields, this term will provide suf- 
ficient restoring force for spin fluctuations to give under- 
damped global spin oscillations. As the magnetic field in- 
creases, the restoring force for such spin waves decreases, 
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eventually vanishing when the synchronized solution be- 
comes unstable. Before this instability occurs, there is 
a transition between underdamped and overdamped spin 
oscillations, and at this transition, the spin wave energies 
have a linear dispersion vs momentum. 

Introducing vortices, the phase locking term favors co- 
alignment of vortices in the two spin polarizations, but 
at sufficiently large A anti-aligned vortex pairs may be- 
come stable. The phase- locking term also makes it possi- 
ble to obtain a stationary complex of rarefaction waves. 
The experimental realization of such solitary vorticity- 
free waves would demonstrate another aspect of super- 
fluid behavior 57,58 in the incoherently pumped polariton 
system. 

Considering the transition between synchronized and 
desynchronized states in the inhomogeneous system, for 
sufficiently small pump spots, the transition is similar to 
the two-mode model. For larger pump spots however, 
an instability to vortex-lattice formation may preempt 
desynchronization, driven by the density imbalance re- 
quired to sustain a synchronized solution. For both large 
and small pumping-spot radius, the phase portraits near 
the critical detuning can be more complicated, show- 
ing synchronized limit cycles with winding number 0, 
desynchronized limit cycles with winding number 2, and 
chaotic behavior. 

In conclusion, the results presented here illustrate the 
wide variety of dynamical behavior that can arise in 
spinor polariton condensates, and suggest that experi- 
mental efforts to investigate such behavior should be fea- 
sible. 
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APPENDIX A: SPIRAL VORTICES IN THE 
SINGLE COMPONENT CONDENSATE 

This appendix will discuss spiral vortices in the one- 
component version of Eq. ([2]), i.e., neglecting J and A. 
It was shown in Ref. [35| that below the critical radius of 
pumping for the formation of a vortex lattice there exist 
states with one or a few spiral vortices. By combining 
vorticity with pumping and decay, one has both radial 
and azimuthal supercurrents, both of which modify the 
density profile, and so these currents interact. This is 
shown in Fig. [20l which shows both the radial velocity 
(main figure) for vortex solutions, and the density profile 
(inset). The vortex solution requires vanishing density at 
the origin, this makes the region at small radii a region of 



u(r) 




0.1 0.2 0.3 0.4 r 



FIG. 20: (Color online) Radial velocity, u(r) — </>'(r), of the 
vortex solutions of Eq. (|A1[) for a — 1, a = 0.01. The numbers 
next to the lines indicate the winding numbers of the vortices 
(1 and 2) and the ground state (0). The TF approximation 
of the gradient flow of the ground state for small a and a 
u(r) = — rcr(/jL — r 2 )/6, is given by the dashed line. The inset 
compares the density, p(r), of a vortex solution for a = 4, a = 
0.27 (fi = 30.28, gray (red) line) with the density of the GPE 
vortex (a = 0,cr = 0) with the same number of particles 
(2tt Jl^l 2 rdr = 1000, fi = 25.68, black line). Both vortex 
solutions have topological charge 1. 



net gain, and leads to an outward- flowing current. Thus, 
in the vortex solutions, there exists both a local maxi- 
mum of radial current, and a point at finite radius (of 
the order of the healing length) at which this current 
vanishes. 

In polar coordinates (r, 0) the wavefunction of a spiral 
vortex of topological charge s in the one-component po- 
lariton condensate takes the form ijj v = f(r) exp[i</>(r) + 
s0], with the equations governing / and u(r) — (j) f (r) 
given by 

( rfu j = (a6(r - r) - <x/ 2 )/ 2 , 

/" + -/'+U-« 2 -^-/ 2 -r 2 / = 0. 
r \ r z J 

Around the center of the trap f(r) and u(r) can be found 
recursively in the form of the power series 

CO CO 

f{r) = ^2a i r^ 2i ~ 1 \ u{r) = J^r 24-1 . (A2) 

2=1 i=l 

To the leading order u(r) ~ a/2(s + l)r showing that 
the stronger the pumping the larger the outward veloc- 
ity is in the vortex core. Only vortices of topological 
charge 1 are dynamically stable. For these, a\ w <5/i, as = 
— ai/i/8..., &2 = (a/i — 8aa 2 )/48..., where the numerical 
value of the vortex core parameter S w 0.583 was first 
calculated by Pitaevskii 59 for the straight line vortex of 
uniform Gross-Pitaevskii condensate. 
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APPENDIX B: IRREGULAR VORTEX LATTICES 

We note that around and above critical detuning, so- 
lutions to Eq. ([2]) may exhibit irregular density profiles 
and irregular vortex lattices. This happens most fre- 
quently when solutions are found starting from Gaussian 
initial conditions, but it may also happen for large tq 
when A is increased stepwise. As an example, we show 
a solution obtained from Gaussian initial conditions at 
J = 1.0,ro = 6.0, A = 6.0 in Fig. EQ This solution 
shows an irregular vortex lattice, reminiscent of turbu- 
lent behavior. 



In many cases, Eq. (|C1|) with s = is sufficient to deter- 
mine ii. However, when there is a vortex at the center of 
the cloud, the numerator and denominator both vanish, 
whereas one of s = ±1 will give a well defined value. In 
such cases, it is then necessary to use Eq. (|C2|) and the 
integral for (L z ) to eliminate Q. 

For plotting the phase portraits, it is necessary also to 
calculate the associated phase, <\> such that = fi. This 
can similarly be defined as: 



Im < In 



i/j(r) d 2 



(C3) 
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FIG. 21: Density of the left- and right-polarized components, 
respectively, showing an irregular vortex lattice. 



APPENDIX C: CHEMICAL POTENTIALS, 
PHASE PORTRAITS AND NUMERICAL 
PRECESSION 

In order to numerically find the chemical potential of 
a condensate with a rotating vortex lattice, it is useful to 
use: 

/i-2*»= J r (CI) 



along with the equation: 



However, when there is a central vortex, the same dif- 
ficulty will arise. To avoid this problem, two possible 
approaches are used. In some most cases it is sufficient 
to calculate = hl — Mi?, and numerically integrate this 
equation to find 0{t). Numerical errors in evaluating hl,r 
can however introduce unphysical precession into the plot 
of 0, 0. To avoid this, in those cases where the nature of 
the portrait is simple, we use 0i+\ = 0^ + x dt x / , with 
f — 1 adjusted to remove the spurious precession. This 
method has been used in Figs. [T2], [13] (right hand panel) 
and [H 

Where the phase portrait involves finer structure, such 
as the middle panel of Fig. [13l the phase portrait has 
instead been obtained by: 



= Im { In 
In 



J e is ^ R {v)d 2 r 
e is ^L(r) d 2 r 



(C4) 



with s = 0, ±1 according to whether there is a central 
vortex. Since the two components are corotating in this 
case, such a definition satisfies: 

= (fji L - 2sQ) - (fjL R - 2sQ) = fi L -m (C5) 

as required. 
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